Winding up by a quench: Insulator to superfluid phase transition in a ring of BECs 
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We study phase transition from the Mott insulator to superfluid in a periodic optical lattice. 
Kibble-Zurek mechanism predicts buildup of winding number through random walk of BEC phases, 
with the step size scaling as a the third root of transition rate. We confirm this and demonstrate 
that this scaling accounts for the net winding number after the transition. 
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Introduction. — In a second order phase transition, 
the critical point is characterized by divergences in the 
correlation length and in the relaxation time. This crit- 
ical slowing down implies that no matter how slowly a 
system is driven through the transition its evolution can- 
not be adiabatic close to the critical point As a result, 
the state after the transition is not perfectly ordered: it 
is a mosaic of domains whose size depends on the rate of 
the transition. This scenario was first described in the 
cosmological setting by Kibble [l[ who appealed to rela- 
tivistic causality to set an upper bound on domain size. 
The dynamical mechanism that determines domain size 
in second order phase transitions was proposed by one 
of us [2| ■ It is based on the universality of critical slow- 
ing down, and predicts that average size of the ordered 
domains £ scales with the transition time tq as £ ~ Tq, 
where w is a combination of critical exponents. This 
Kibble-Zurek mechanism (KZM) for second order ther- 
modynamic phase transitions was confirmed by numer- 
ical simulations Q and tested by experiments in liquid 
crystals N| , superfluid helium 3 [5| , both high-T c [6| and 
low-T c [7f superconductors, and even in non-equilibrium 
systems [8j . With the exception of superfluid 4 He - where 
the situation remains unclear Q , experimental results are 
consistent with KZM (see [l(| for a review) . Spontaneous 
appearance of vorticity during Bose-Einstein condensa- 
tion driven by evaporative cooling was recently reported 
11 1 . This confirms KZM predictions [12 |. and is further 
elucidated by numerical studies of BEC formation [l]| • 

Our goal is to study dynamics of a quantum phase 
transition in a simple yet non-trivial example that can 
be implemented experimentally. Quantum phase transi- 
tions we consider differ qualitatively from finite tempera- 
ture transitions. Most importantly, evolution is unitary, 
so there is no damping, and no thermal fluctuations to 
initiate symmetry breaking. Recent work on the dynam- 
ics of quantum phase transitions is mostly theoretical, 
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possible exception: Ref. [24j on the transition in a spin-1 



BEC. Generic outcome of that experiment is a mosaic of 
ferromagnetic domains whose origin was attributed to a 
sudden quench limit of KZM. This explanation is sup- 



ported by theory 25 1. 

Model. — Bose-Hubbard model is a paradigmatic 
example of a non-integrable quantum critical system. It 
describes cold bosonic atoms in an optical lattice [2^] . In 
dimensionless variables, its Hamiltonian reads 



N N 
s=l ~ 



(1) 



Here N is the number of lattice sites and n is an aver- 
age number of atoms per site. This model with periodic 
boundary conditions (which we assume) should be di- 
rectly experimentally accessible in a ring-shaped optical 
lattice [301 ] ■ For an integer n, the transition from the 
Mott insulator (small J) to the superfluid phase (large 
J) is located at J c ~ n~ 2 [2&| . 

We drive the system through its critical point by a 
linear quench with a quench timescale tq: 



J{t) = t/TQ 



(2) 



In an experiment one can increase Josephson coupling J 
by turning off the optical lattice potential as in [29( . The 
initial state is the Mott insulator ground state at J = 0, 
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with the same atom number at each site. We assume n 3> 
1: This large density limit is accessible experimentally. 

Numerical approach. — We replace annihilation 
operators a s by complex field <j> s , a s w \fn <f> s , which 
is normalized, J2 s =i l^s| 2 = N , and evolves with the 
time-dependent Gross-Pitaevskii equation 

i—— 

dt 

These approximations are accurate for n — > oo, when the 
critical point J c ~ n~ 2 — > 0. 

In the truncated Wigner method we employ quantum 
expectation values are given by the averages over stochas- 
tic realisations of the field (j) s (t) [H HI, 123] ■ For example, 
the correlation function becomes 
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Here (..) means quantum expectation value while the 
overline is an average over realizations. All realizations of 
<j> a (t) evolve with the same deterministic Gross-Pitaevskii 
equation ([¥]), but they start from different random initial 
conditions which come from a probability distribution 
depending on an initial quantum state. The initial Mott 
state ([3]) corresponds to initial fields 
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with independent random phases 6 S <= [0, 2%): The Mott 
state has the same number of n particles at each site 
(i.e., 1 0s (0)| = 1), and, hence, indeterminate phases, that 
translate into random 8 S . 

Kibble-Zurek mechanism. — In an optical lattice 
with BEC pools that become gradually connected with 
Josephson couplings in accord with Eq. ([2]) it is natural to 
rephrase KZM: Rather than seek distance £ over which 
phase remains more or less the same we compute size 
A9 S of a typical phase step between neighboring sites. 
One could use it to deduce the size of domains £ over 
which winding number changes by one, and get the ac- 
cumulated phase from square root of circumference of the 
whole ring of BEC pools measured in units of £, as in Q • 
However, the same result obtains from a random walk 
between neighboring sites, with the corresponding step 
size A8 S . We now compute A9 S as a function of tq. 

The Gross-Pitaevski equation (fj| can be lin- 
earized in small fluctuations 8<j) s around uniform 
large background, <fi s — 1 + 5<ft s , and S(f> s can 
be expanded in Bogoliubov modes as 5<ft s — 
J2k (bkUke lks + frfc^e~ lfcs ) with pseudomomentum k. 
For constant J we have bk(t) = bk(0)e~ luJkt with cok = 
2-J J(l — cos fc) [1 + J(l — cos kj] and stationary Bogoli- 
ubov modes u k = —Nk [1 + 2J(1 — cos A:) + Wk] , Vk = 



Afk where Afk are such that it? - 



1. In the Josephson 



regime, when J -C 1, we have Vk ~ — Uk so that purely 
imaginary 6<p s in 4>s = 1 + S<p s is a phase fluctuation. 
However, for our random initial conditions ©, this lin- 
earization is justified only for short wavelength modes of 
<j) s , with k ~ ±7r, for whom the modes with longer wave- 
length are a locally uniform large background. From now 
on we focus on the short wavelength modes because they 
determine variance of the nearest-neighbor A9 S . 

When k w ±ir and J« 1, then uo k w 2\f2J . Early 
in the linear quench ([2]) this lo^ is small, so that early 
evolution of the short wavelength modes is approximately 
impulse i.e. their magnitude remains the same as in the 
initial Mott state and, consequently, A# 2 ~ 1 in this 
impulse stage. The impulse approximation breaks down 
at J ((3) when the transition rate uik/uik equals Wfe, 

uj k /uk ^ ^fe , (7) 
and evolution becomes adiabatic. Eq. Q leads to 
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(8) 



which is consistent with J< 1 when tq ^> 1. 

The crossover from impulse to adiabatic evolution at 
J is the key ingredient of KZM. In the following adia- 
batic evolution after J but before J « 1, short wave- 
length phase fluctuations scale as 54> s ~ J^ 1 / 4 because 
the mode amplitudes \b k \ do not change, but u k and 
Vk follow stationary Bogoliubov modes u k sa — v k ~ 
— 1/2(2J) 1 ^ 4 . Consequently, A9 S has variance scaling as 

at J that A6*J ~ 1, phase fluctuations must shrink as 

.-1/3 



Given the boundary condition 
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J -1 / 2 while J« 1. 
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On the other hand, when J> 1 then stationary modes 
Uk ~ 1 and v k ~ do not depend on J and A# 2 does not 
depend on J either. This means that A# 2 must stabilize 
between the regimes of J -C 1 and J ^> 1 i.e. around 
J ~ 1 where it takes its final value 
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(9) 



which scales with a power of w = 1/3. 

This variance determines e.g. the correlator C\ in 



1 - d = 1 - cos A6» s 
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(10) 



for tq > 1. Kinetic hopping energy per particle K\ 
is expected to stabilize for J ^ 1, when the hopping 
term dominates over the non-linearity in Eq. (U) and Ki 
becomes an approximate constant of motion, see Fig. [TJ 

Key ingredients of KZM are confirmed by our simula- 
tions: Phase performs a random walk that is markovian 
to a good approximation. Moreover - as seen in Fig. [T]- 
its size is consistent with the above predictions. 

Winding number. — Condensate wavefunction is 
single-valued. Therefore, phase accumulated 0# after 
R = N steps defines integer winding number: 



N 



(ii) 



where Arg(...) 6 (— it, 7t]. A random walk of phase, with 
the variance of nearest neighbor phase differences scaling 
as in Eq. ([§]), gives winding numbers with variance 



N T. 
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Q 



(12) 



There are two limits where this scaling is bound to fail. 
For very fast quenches with tq <C 1 phases are completely 
random between neighboring sites, so A6* 2 = 7r 2 /3, and 
W$f = N/12. For quenches so slow that Wjj < 1 the 
na ture of the problem changes, leading to steeper falloff 
of W% with tq 0,13. Between these two limits the i- 
scaling in Eq. (112|) for the winding number is confirmed 
by our numerical results in Fig. [21 
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FIG. 1: Kinetic hopping energy K\ = 1 — Ci » 1 — cos A6 S ps 
^A^s as a function of rescaled J/ J for different tq is seen in 
A. When J <C 1, all the plots overlap demonstrating that J = 
i~q 2 ^ 3 is the relevant scale for J C 1. Individual plots depart 
from this small-J "common bundle" at J ~ 1, or J/tq 2 — 

2/3 

Tq , when Ki = 1 — Ci is expected to stabilize. In B, we show 
Kr = 1 — Cr at J — 10 as a function of tq for R — 1, 5. 
Data points for each J? were fitted with lines, their slopes 
giving exponents close to the |-scaling predicted in Eq. (|10[) 
with error bars on their last digits. Size of typical phase step 
can be estimated as AOs — Kr/R when Kr <C 1, and already 
this rough estimate yields a good approximation of winding 
numbers shown in Fig. [2] 

Correlation function. — Constant amplitude and 
Gaussian distribution of phase Qr after R steps imply 

where or is dispersion of Qr = X)fLi Arg {4>s+i<frt) 
which after R = N steps becomes the winding num- 
ber in Eq. lfTTj) i.e. Wn = Qn/^tt. For a random walk 
a R = RA9^, which leads one to expect: 

C R ~ exp (-AA0J/2) = exp(-i?/0 , (14) 

Using Eq. ((9} we would expect scaling £ ~ Tq 3 '. 

Numerical simulations confirm exponential correla- 
tions, see Fig. [31 but correlation lengths £ measured at 
J = 10 are better fitted by £ ~ t q 45 - O n the other hand, 
early on in the quench, for smaller values of J <C 1, 
correlation length exhibits £ ~ t q^ 3 - ^ seems that in- 
termediate scales are subject to phase ordering between 

2 /3 

the freezeout at J ~ Tq and the final J = 10. Sim- 
ilar post-transition phase ordering was observed in the 
integrable quantum Ising chain [2 ll | . 

On the other hand, winding number continues to scale 
with Tq , see Fig. [2] It is not too surprising that it is 
insensitive to phase ordering: While in our simulations 
winding number is not really stable following the freeze- 
out, it changes much less frequently than smaller scale 
excitations, as its topological nature leads one to expect. 
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FIG. 2: Variance of winding number measured at J = 10 
as a function of tq for lattice sizes N — 512, 256, 128. Here 
point sizes equal error bars. The data points with tq > 2 and 
> 2 were fitted with the solid lines giving slopes close 
to the predicted — | in Eq. (|12|l . is shown over a wider 
range of tq to show the saturation for nearly instantaneous 
quenches, when tq < 2, and the crossover to steeper slope 
when Wfj < 2. In the inset we show a rescaled W^/N for 
N = 512, 256, 128, 32, 8 to demonstrate that W% ~ N in the 
KZM regime of tq > 2 and W% > 2. 

Summary. — We have investigated the process mak- 
ing a single condensate wavefunction out of many — N - 
independent BEC pools. We conclude that, in the ring 
geometry, the overall winding number Wn (which will 
set up persistent current) can be predicted using simple 
idea of a random walk in phase between the initially in- 
dependent BEC fragments For very quick quenches 
this leads to saturation at Wn = N/12. Slower quenches 
lead to scaling of Wn with the rate of reconnection that 
can be inferred from the Kibble-Zurek mechanism. 

Correlation functions also exhibit behavior consistent 
with a random walk in phase. Initially, correlations scale 
in a way that is directly related to healing length at the 
instant when dynamics of the system becomes faster than 
the rate of change of its Hamiltonian However, while 
winding number "remembers" this scaling as Joseph- 
son couplings increase, correlations on smaller scales 
evolve. In thermodynamic transitions similar phase or- 
dering associated with diffusion is responsible for the 
post-transition smoothing of the order parameter struc- 
ture, so that - eventually - only topological defects still 
"remember" initial state of the system. In our model 
evolution is completely reversible. Therefore, diffusion 
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FIG. 3: Correlation function Cr for different tq measured 
at J = 10 on a lattice of N = 512 sites. This logarithmic plot 
demonstrates that Cr is exponential in R. 



cannot smooth out small scale structures. However, evo- 
lution itself appears to redistribute energy between the 
excitations. This may be regarded as a quantum ana- 
logue of phase ordering. Correlations on intermediate 
scales change, but (as was also the case in thermody- 
namic phase transitions) small-scale evolution does not 
affect the topologically protected winding number Wn- 

Our model ignores decoherence and damping that are 
likely to intervene in the laboratory experiments with, 
say, gaseous BECs. It is relatively easy to modify equa- 
tions and introduce damping "by hand" . There is how- 
ever no unique prescription for it (although one could 
appeal to presence of a dilute thermal cloud, as in simple 
models of BEC decoherence [3l|). In experiments dissi- 
pation and decoherence are inevitable. We expect dissi- 
pation to affect small scales, but leave the topologically 
conserved Wm intact. This is based on a limited number 
of simulations we have conducted where different models 
of dissipation were tried out. Above all, this is corrobo- 
rated by the experiment 11] where sudden reconnection 
of N — 3 uncorrelated condensates led to relaxation to a 
condensate with stable vortices - stable winding number. 



It is also consistent with the recent numerical results [32 1 . 
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